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Abstract 

This paper investigates the dynamics of biomass in a marine ecosystem. A stochas- 
tic process is defined in which organisms undergo jumps in body size as they catch 
and eat smaller organisms. Using a systematic expansion of the master equation, 
we derive a deterministic equation for the macroscopic dynamics, which we call the 
deterministic jump-growth equation, and a linear Fokker-Planck equation for the 
stochastic fluctuations. The McKendrick-von Foerster equation, used in previous 
studies, is shown to be a first-order approximation, appropriate in equilibrium sys- 
tems where predators are much larger than their prey. The model has a power-law 
steady state consistent with the approximate constancy of mass density in logarithmic 
intervals of body mass often observed in marine ecosystems. The behaviours of the 
stochastic process, the deterministic jump- growth equation and the McKendrick-von 
Foerster equation are compared using numerical methods. The numerical analysis 
shows two classes of attractors: steady states and travelling waves. 
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1 Introduction 



Marine and freshwater ecosystems exhibit a remarkable regularity in the relation between 
abundance of organisms and their body masses. Treating organisms simply as particles 
of different size, i.e. ignoring taxonomic identity, the total biomass (abundance x body 
mass) in logar i thmic intervals of body mass is o bsery e d to be approximately constant 



(jSheldon et al. 



1972 



19771 ; iBoudreau and Dickid . Il992j ; Kerr and Dickie! 12001 ). Equiva- 



lently, the logarithm of abundance expressed as a function of the logarithm of body mass, 
often referred to as a size spectrum, is approximately linear with a gradient close to — 1. 
Removing the logarithms, this is equivalent to density in mass space being a power func- 
tion of mass with an exponent —2. This approximate regularity applies over a wide range 
of body size from micro-organisms to large verteb r ates, and has been the subject of much 



research and discussion in ecology (j Sheldon et al. , 1972 ; Piatt and Penman . 1978 ; Heath 



1995 



Marquet et all 12005( 1 . 



Understanding of the dynamics of biomass flow that lead to this regularity is important: 
the biomass of most marine ecosystems supports major fisheries that play a significant role 
in the economies of maritime countries. The dynamics are often studied by means of a 
partial differential equation (PDE), in which abundance is taken as a fun ction of both body 
mass and time. The PDE is motivated by a model of McKendrick (119261 ) and von Foerster 
(119591 ) . in which abundance is a function of age and time. We will follow the convention 
of calling this PDE the McKendrick-von Foerster equation, notwithstanding the change in 
variable from age to size. 



The McKendrick-von Foerster equation was first adopted by Silvert and Piatt ( 119781 ) 
in a model allowi ng gr owth and mortality to be functions of body mass. Following this, 
Silvert and Piatt (119801 ) coupled growth at one size to death at another, because organisms 
grow in size spectra by eating smaller organisms. More recently, the approach has been 
exten ded, first to allow organisms to eat those at all smaller sizes (jCamacho and Sola 
200ll ). and second, by using a feeding-kernel fun ction, to allow them to eat organisms in 
a restricted size range (IBenoit and Rochetl . 120041 ) . PDEs of th is kind are now being used 
quite extensively to understand proc e sses in marine ecosystems ( Andersen and Beyer . 20061 : 



Maury et al. 



2007 



Andersen et all . 120081 ). It can, for instance, be shown in numerical 



analyses that the PDE at stea dy state gives size spect ra with gradients which are similar 
to those in marine ecosystems (IBlanchard et all 120091 ) . 



The McKendrick-von Foerster equation is implicitly assumed to be an appropriate ap- 
proximation for an underlying stochastic process in which individual organisms grow by 
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eating prey items. A first investig ation of the relationship between the PDE and the 



stochastic process (ILaw et all 120091 ) showed that the PDE could describe the approach to 
a steady-state size spectrum. However, the stochastic process could also develop travelling- 
waves; although these were also possible in the PDE, the properties of these waves were 
somewhat different. The research described in the present paper was motivated by these 
discrepancies. 

A possible source of these discrepancies is that the McKendrick-von Foerster equation 
was originally conceived of as a model for organisms indexed by age, rather than by weight. 
Age and weight do not change in quite the same way over time. An organism grows older 
continuously, whereas its weight grows in jumps each time it finds a prey item to feed 
upon. As time progresses, organisms which start at the same age clearly remain the same 
age as each other, whereas organisms which start at the s ame weight in general do not 



remain the same weight as each other. Pfister and Stevens (120021 ) stressed the impo rtance 
of growth variability in cohorts of organisms. Motivated by this, Gurney and Veitch ( 120071 ) 
considered the consequences of allowing growth to be both a random variable and also size- 
dependent, in a von Bertalanffy growth model. However, the emphasis in dynamic size 
spectra is somewhat different because variation in body weight here emerges from random 
encounters with prey items of various weights. 

In this paper we therefore start from a stochastic process in which organisms undergo 
jumps in body size as they catch and eat smaller organisms. We introduce this individual- 
based stochastic process in 12.11 and describe it as a population-leve l model in Section 1 2.21 
In Section [2~51 we use a systematic expansion of the master equation (Ivan Kampenl . ll992l ) to 
derive an equation for the macroscopic dynamics (which we call the the deterministic jump- 
growth equation (fT2"j) ) and a Fokker-Planck equation for the stochastic fluctuations. We 
also provi de an appendix with an alternative derivation of a Langevin equation, following 



Gillespie (120001 ). to clarify an issue unresolved by the systematic expansion. Section 1231 
shows that the McKendrick-von Foerster equation is a first-order approximation of the 
deterministic jump-growth equation, which applies at steady state when predators are 
much larger than their prey. In Section 12.61 we show that our model has a power-law 
steady state and we derive an approximate analytic expression for its exponent, thereby 
showing that the steady state is consistent with the approximate regularity seen in marine 
ecosystems. However, the steady state is not necessarily an attractor. In Section [3] the 
behaviour of the deterministic models and of the stochastic model are com pared using 



nume rical methods. As in the case of the McKendrick-von Foerster equation (ILaw et al 



20091 ). certain parts of parameter space allow a travelling- wave solution. 
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2 A dynamical model of size-dependent predation 



2.1 An individual-based stochastic process 

We model predation as a Markov process. The primary stochastic event comprises a 
predator of weight w a consuming a prey of weight w& and, as a result, increasing to become 
weight w c (Figured]). Predation is inefficient and, in keeping with ecological convention, 
we assume that a fixed proportion K of prey mass is assimilated by the predator so that 
u! c = w a + Kwb (the assumption of constant K could be relaxed). We call this model the 
'jump-growth model' because the changes in the weight distribution are caused by discrete 
steps in body size as predators eat prey, and the mortality that comes with this predation. 

It would be easy to add additional events to the jump-growth model to account for 
natural death and for birth (recruitment) but, as we will see, for the purpose of this paper 
of explaining the observed power law size spectrum these additional events are not required, 
and we will therefore restrict our attention to the pure predation events. 

The next three subsections will be concerned with the derivation of equations describing 
the time evolution of the weight distribution that follows from this stochastic process. The 
main result from these sections that we will use further in this paper is the deterministic 
jump-growth equation (fl2|) given in Section EH That equation has an intuitive explanation 
in figure [T] and is enough to follow the remainder of the paper. 

A mathematically rigorous treatment of the individual-based model may be possible fol- 



lowing th e techniques deve 



example ( IFinkelshtein et al 



oped for stochastic processes on configuration-space, see for 



20091 ). In this paper we will instead pursue a heuristic treat- 



ment of a corresponding population-level model. 



2.2 A population- level master equation 

Instead of keeping track of the weight of each individual, we aggregate individuals of similar 
weights into weight brackets, and follow the number of individuals in each bracket. We 
introduce a set of weights Wi and corresponding weight brackets [wi,Wi + i), with i £ Z. In 
practice, the size of the weight brackets should be chosen small enough so that discretisation 
errors are small. The weight distribution of organisms in a large fixed volume Q is described 
by a sequence of numbers "[..., ri-i, uq, ni, . . . ], where Hi is the number of organisms in Q 
with weights in the i-th weight bracket between Wi and w i+ i. 
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Figure 1: The primary predation event replaces an individual predator and prey by a new, 
larger predator individual. Taking some arbitrary weight w, there are two ways in which 
an individual can change from this weight: by feeding and thereby increasing in weight, 
and by being eaten and so disappearing altogether. There is also one way in which an 
individual can become weight w. by being of smaller weight and feeding on a prey of just 
the right size to become weight w. These events are reflected in the three terms of the 
deterministic jump-growth equation f[T2l in Subsection 12.41 



Let kij/Q denote the rate constants for the predation events, where the indices of k 
are ordered: predator before feeding, prey. Then the probability in an infinitesimal time 
interval dt for any one of n« organisms in weight bracket i to eat any one of the rij organisms 
of weight bracket j is kijVL^niUjdt. The time evolution of the probability P(n, t) that the 
system is in the state n at time t is then given by the master equation 

9P( £' t] = E if ^ + vfa + x ) p ( n - u *> *) - *)] . (!) 

where n — = (. . . , rij + 1, . . . , rii + 1, . . . , n\ — 1, . . . ), and / is the index of the weight 
bracket wi < Wi + Kwj < wi + \. The first (positive) term in flTJ corresponds to having one 
extra predator in weight bracket (i), one extra prey in (j), and one less predator in (/), so 
that one predation event will move the system from state n — Uij into state n. The second 
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(negative) term corresponds to another such predation event that moves the system out of 
state n. Hence the master equation is commonly referred to as a "gain-loss" equation. 



2.3 Separation of macroscopic behaviour and fluctuations 

The master equation ([1]) has non-linear coefficients and is difficult to solve analytically. 
We therefore make use of the property that, in systems of sufficiently large volume Q, the 
fluctuations are relatively small because they are suppressed by a factor of the square root 
of Q; the conditions required for this to be tr ue are given in Appendix HI In this section 



we adopt the approach of van Kampen (119921 ) . carrying out an expansion of (CQ) in Q, and 



collecting together the highest-order terms in Q. To do this, it helps to rewrite the master 
equation ([T]) using a step-operator notation: 

= E k ~i (wsr 1 - l ) ( *)) • ( 2 ) 

Here E, is a step operator that acts on a function f(n) to give Ej/([. . . , rij, . . . ]) = 
/([..., 7ij + 1, . . . ]); similarly Ej acts on a function f(n) to give Ej/([. . . , rij, . . . ]) = 
/([..., rij + 1, . . . ]); conversely E 2 _1 acts on a function f(n) to give E / ~ 1 /([, . . . , n h . . . ]) 
= /([..., rii — 1, . . . ]). Thus (J2J) is just an alternative notati on for ([!]). For further expla- 



nation of the step-operator notation, see van Kampen (119921 page 139). 



Following the method used by van Kampen ( 119921 ). we separate each random variable rii 
into a deterministic component <fii(t) which describes the density of individuals in weight 
bracket i, and a random fluctuation component as 

ni = fl&(*) + fi*&(t). (3) 

On average the number of individuals will be proportional to the system size Q, by the law 
of large numbers, and that is the reason for the factor of Q multiplying <pi(t). That the 
fluctuations are proportional to the square root of the system size should be justified by 
some sort of central limit theorem. A heuristic justification is given in appendix HI Thus 
disaggregating rii in this way leaves two variables 0« and £j which no longer scale with 
the system size. We assume that Q is so large that the discrete nature of n is no longer 
noticeable at the level of <f> and $ and we can treat them as continuous variables. 

The new random variables are described by a probability distribution n(£, £) = 
r2 1 / 2 P(n,t). An equation for the time evolution of this probability distribution is ob- 
tained by substituting the change of variables ([3} into the master equation (J2J). Care needs 
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to be taken because this change of variables is time-dependent. This has the consequence 
that 

dP(n, t) _ 1/2 dnjj, t) v dfr dU(t t) 
at ot 2-j dt d& ' u 

Here we used the property that Vt~ 1 ^ 2 d$ t / dt = —d<f>/dt when we keep n fixed. The operators 
Ej which change to + 1 now change to £j + fi" 1 / 2 and can therefore be expanded as 



Substituting all this into the master equation (J2J) gives an equation with terms containing 
various different powers of the system size Q. 

The highest order terms are at order Q°. They only contain the macroscopic variables 
(pi and vanish if these satisfy the deterministic equation 

= ^2 (.-kij<Pi<Pj - kji^j^i + k mj (j) m (pj) , (6) 

3 

where m is an index for the weight bracket: w m < W{ — Kwj < w m +i- The three terms in 
(EJ) are in keeping with the intuition given by Figure [TJ Losses from weight bracket i (the 
negative terms) occur because individuals in this bracket eat prey and become heavier, and 
because these individuals are themselves eaten. Gains into weight bracket % (the positive 
term) occur through smaller predators growing into this bracket by eating prey. Imposing 
the deterministic equation fl6]) is not the only possible way to make the terms of order Q° 
vanish, but it is the most natural and is independently derived in appendix HI 

Terms at the next order, f2~ 1//2 , give the linear Fokker-Planck equation for the probability 
distribution II (£) of the fluctuations, 

dU \ - , d „. 1 x - „ d 2 „ 

« =-E^te n ) + 2E B «ap|- n . P> 

where the coefficients and are independent of the fluctuations If we introduce 
the objects fc^-j and fak by 



k; 



kij if wi < Wi + Kwj < wi + i 



\ otherwise 

fill = g (^Ji + M ( 9 ) 
then we can give the succinct expressions 

An = ^2 fwfai Ai i = (fwfa ~ /«t0/) > ( 10 ) 
ji i 

B-n = ^2 fju&j&h B ij = ^ (fyifa&j ~ fujfcfa - finMi) ■ (H) 
jl I 
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Because the fluctuations are damped by a factor of f2 1//2 , in the remainder of this paper 
we concentrate on studying the deterministic equation fl§]). 



2.4 The deterministic jump-growth equation 

For analytical calculations and also for conceptual considerations it is convenient to work 
with the continuum limit of the macroscopic equations ([6]). This limit is obtained by 
writing the size of the weight brackets as Aj = u>i+i — Wi = //jA and taking the limit 
A — > 0. The discrete set of variables <pi is replaced by a continuous density function (j)(w) 
satisfying <p(wi) = 0j/Aj. This function <f>(w) describes the density per unit mass per unit 
volume as a function of mass w at time t; it therefore has dimensions M -1 L -3 . The sum 
over weights in fl6]) is replaced by an integral, J2i A* — > f dw. The rate constants kij are 
replaced by a feeding rate k(w,w') satisfying k(wi,Wj) = k^. The macroscopic equation 
(El) becomes 



d(f)(w) 
dt 



— k(w,w')(f)(w)(f)(w') 

— k(w' ,w)<p(w')<p(w) 

+ k{w - Kw', w')(j){w - Kw')(f)(w'))dw' . (12) 

We call this equation the 'deterministic jump-growth' equation. The three terms in (ITS]) 
are equivalent to those in ([6]), and correspond to the idea in Figure (JTJ that there are two 
ways to leave weight w and one way to enter it. The terms represent, in order: feeding 
on prey to become larger than weight w, being fed upon and removed from the system 
entirely, and feeding on prey of exactly the right size to become weight w. 



Following Benoit and Rochet (120041 ) we assume that the feeding rate takes the form 



k(w,w') = Aw a s (w/w') . (13) 

This states that the rate at which a particular predator of weight w eats a particular prey 
of weight w' is a product of the volume searched per unit time and a dimensionless feeding 
preference function s. The volume searched per unit time Aw a depends on the predator's 
body weight w, raised to the constant power a. A is a constant volume searched per unit 
time per unit mass a . The feeding preference function s depends only on the ratio w/w' 
between predator weight and prey weight. In practice this feeding preference function will 
be peaked around a preferred predator:prey weight ratio B. 

When the parameter K, that describes which proportion of the prey mass is assimilated 
by the predator, is equal to 1, the deterministic jump-growth equation (TT2|) reduces to 
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the Smoluchowski coagulation equation ( ISmoluchowskil . Il916l ). that is used to describe 
the clumping together of particles, for example in aerosols. However the rate kernels 
used to describe coagulation differ greatly from our localised feeding rate kernel ffTB"]) . 
Typical choices in the coagulation equation are k(x, y) = x + y or xy or other homogeneous 
expressions and these lead to ver y different beh aviour such as the formation of one giant 
cluster (gelation); see for example lAldousI ( 119991 ). 



2.5 Relation to the McKendrick— von Foerster equation 



The deterministic jump growth equation (j!2p is not the same as the McKendrick-von 



Foerster equation which has 



fSil vert and Plat 



2009 



Law ct al. 



1978 



1980 



jeen widely used to describe the dynamics of size spectra 



Benoit and Rochet 



2004 



Maury et al. 



2007 



Blanchard et al 



20091 ) and which reads 



where D is the per capita death rate at weight w from predation by larger organisms, 



D(w) = I k(w' ,w)(j)(w')dw' , 



and G is the growth rate at weight w from feeding on smaller organisms, 



G(w) = / Kw'k(w,w')(j)(w')dw' 



(15) 



(16) 



Here we show that ffl4l) emerges as an approximation to ffl2|) in the case where the typical 
prey is small in size compared with the predator. Such an assumption is reasonable in many 



cases, because predators tend to be of an order 



prey (I Cohen et al 



1993 



2 to 10 3 times the body mass of their 



Jennings and Mackinsonl . 120031 ) . So the feeding kernel is strongly 



peaked around w' = w/B with B large. Taki ng into account furthe r the i nefficiency with 
which prey mass is assimilated (K fh 10 _1 ) (IPaloheimo and Dickid . Il966l ). there is some 
justification for treating Kw' « w in the last term of (IT2~1) . This motivates a Taylor 
expansion of this term around w, 



k(w — Kw',w')(j)(w — Kw') = k(w,w')(j)(w) 

d 



+ (-Kw')— (k(w,w')(j)(w)) 
aw 

(-Kw') 2 d 2 n , AX 
+ 2! g^( k (w,w'Mw)) + 



(17) 



S 



Substituting this back into (|T2|) gives 




2dw 2 
+ R, 



k{w' , w)4>{w)4>{w')dw' 
- I Kw'k(w,w')(j)(w)(j)(w')dw' 



(Kw') 2 k(w, w')(f)(w)(j)(w')dw' 



(18) 



where the remainder term R is given by 



oo 



(-1)" & 



1 



n\ dw 



(Kw') n k(w, w')(f)(w)(j)(w')dw' . 



(19) 



n=3 



The first two terms in (TIB"]) correspond to those in the McKendrick-von Foerster equation 
(TH|) . For ecosystems near to steady state, where <fi(w) is close to a power law (as we will 
see in the next section), the higher order terms are suppressed by factors of K/B and are 
therefore small. Thus the McKendrick-von Foerster equation is a good approximation for 
f fl2|) near the steady state and when prey are typically much smaller than their predators. 
But the higher order terms are not necessarily small in non-equilibrium ecosystems. In 
particular, the McKendrick-von Foerster equation is a less good approximation if there is 
a travelling wave attractor, see Section 13.21 

One way to understand the difference between (fT2l) and ffT^|) is that f JT2|) models the 
discrete, variously-sized jumps in predator mass as predators feed and grow. This captures 
the property of the stochastic model that individuals, starting at a given weight, are able 
to develop a range of weights over the course of time. In contrast to this, the McKendrick- 
von Foerster equation (fT^j) assumes smooth growth along the weight axis. Spreading of 
body size can be incorporated in (f!4|) by introducing the diffusion term, the third term in 
(ITS]) . The source of such diffusion is the deterministic jump-growth equation (i.e. terms 
of order f2°), so diffusion is attributable to the deterministic, as opposed to the stochastic, 
component of the full process. 

2.6 Steady-state solution 

In marine ecosystems, abundance of organisms within body mass classes averaged over 
space and seasons often changes rather little, suggesting that they may be close to a 
steady state. In such circumstances and when abundance and mass are both expressed 
as logarithms (i.e. as a power spectrum) the relationship is approximately linear with a 
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gradient often close to -1, which implies a power law with an exponent -2 in the untrans- 
formed variables. This leads to the important regularity of marine ecosystems that the 
total biomass is approximately constant when expressed in logarithmic intervals of body 
mass. 



Benoit and Rochet (120041 ) found that the McKendrick-von Fo erster equation has s teady 



state solutions of t he po wer law form 



Camacho and Sole 



w 



oc w 7 (see also 



Piatt and Penman. 1 1978 



20011 ). and we will show that the same is true for the deterministic 



jump-growth equation (1121) . Of course in the real world such a power law will have to 
break down for very small weights (where otherwise the power law would predict an infinite 
number of very small individuals) and for very large weights (where the power law would 
predict a non-zero density of arbitrarily large individuals). Indeed, in a real system with 
a finite number of individuals, a model just having predation events could not have a 
non-trivial steady state because the number of individuals would continue to decrease. 
A non-zero steady state is possible only if there is an inexhaustible reservoir of small 
individuals. In our model the power law spectrum provides this reservoir automatically. 
In a more realistic model one would need to model the plankton as well as recruitment. 

A steady state solution 4>(w) of f[T2"j) has to satisfy the equation 



k(w, w')<j){w)4>{w')dw' 
k(w', w)4>{w')<f){w)dw' 
k{w — Kw', w')(j)(w — Kw')(f>(w')dw' , 



(20) 



+ 



If we substitute the power law Ansatz <f>(w) oc w~ 
for the feeding rate, change to a new integration variable r = w pre dator/wp rey and cancel 
some overall factors, we obtain an equation for the exponent 7, 



into this equation, use the form ([13 



= fin) 




_l_ r 



r + K) 



-a+27-2 



dr. 



(21) 



The existence of a power law steady state can now be proven using the same argument 
as that given by Benoit and Rochet (120041 ) in the case of the McKendrick-von Foerster 
equation^. The argument goes as follows. If we assume that predators are bigger than their 
prey, then for 7 < 1 + a/2, f(j) is less than zero. Also, f(j) increases monotonically for 
7 > 1 + a/2, and is positive for large positive 7. Therefore there will always be a unique 



7 for which f(j) is zero and thus a unique steady state of the form 



lw) oc w 



-t. If we 



2 We thank one of the referees for pointing this out. 
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allow predators to be smaller than their prey, situations with no power law steady state or 
multiple power law steady states can be found. 



The numerical value of the power law exponent 7 is of particular interest because 7 is 
known to have a value close to 2 in marine ecosystems (see Section [1]). In the special case 
that K — 1, and a = 1, a value 7 = 2 does in fact satisfy f[2~Tj) . A value of a close to 1 is 
biologically reasonable as this means that the volume searched by a predator is proportional 
to its body weight (see Equation (1X51)). a l thoug h the limited information available suggests 



a value slightly lower than a = 1 ( Ware . [1978). More generally, 7 = (3 + a)/2 will satisfy 
( l2"Tj) for any a, with K = 1. 



A value of K close to 1 is unrealistic: K m 0.1 would be more appropriate (IPaloheimo and Dickid . 



19661 ) because only a small proportion of food ingested is assimilated into extra body weight. 
To treat this case analytically we make the assumption that predators feed only on prey of 
their preferred size, i.e., we set the feeding preference function in (121]) to the delta function 




s(r) = S(r — B). In that case ( |2T|) reduces to 

= -B^- 2 - B a -^ + B a ^(B + K)- a+2 ^ 2 . (22) 
A Taylor expansion in K/B gives 

0^(2 7 -a-2)-^-fi- 27+Q+2 , (23) 
B 

and the Lambert W function can be used to express 7 explicitly as a function of the other 
variables 

1 / ur (R 1™ \ 

(24) 

At K = 1 and a = 1, (EH produces 7 = 2 because W(BlogB) = log B. For K < 1, 
the exponent 7 increases as either K or B decrease, because in either case less mass is 
transferred to larger organisms. Notice however that the dependence of 7 on K and B is 
weak; for instance, if K = 0.1 and B = 100 (still with a = 1), the exponent only increases 
to 7 = 2.21. Thus if K and B are given biologically reasonable values the steady-state of 
the model is broadly consistent with the empirical property of marine ecosystems that 7 
is close to 2. 

The ecological literature contains a relationship between the parameter 7, and K and B 
based on a quite different premise, that the metabolic rate of organisms scales with body 
weight as w 3 ^ 4 . It can be shown from this scaling that 

3 log if 

7 = 1 + — 25 

' 4 bgB v ; 
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in the absence of any consideration of dynamics (IBrown et ajj , |2QQ4j ). There is some resem- 
blance between this equation and (J23J), which becomes evident from taking the asymptotic 
approximation for the Lambert W function 



W(z) = log z — log log z H 



(26) 



in (J21J), giving an expansion in which the leading terms are 
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log-K" _ log log (f log B) log log B 
log B log B log B 



+ 



(27) 



Both fl25|) and (1271) contain the term (log K)/ (log B), but are not the same. From a 
biological standpoint the equations have the important difference that the relationship in 
( )2~Tj) follows simply from dynamical bookkeeping o f biom ass, without any assumption about 



We stress that, although some properties of the steady state have been described here, 
we have not investigated analytically the region of parameter space in which the steady 
state is an attractor. The next Section shows by means of numerical methods two 
classes of attractor: a steady state of the kind described above and a travelling wave. 

3 Numerical results 

Here we use numerical methods to compare some properties of the stochastic jump-growth 
model 02]), the deterministic jump-growth equation ( fl2|) and the McKendrick-von Foerster 
equation ( IT3|) . 

Body sizes can span at least ten orders of magnitude in real ecosystems, and it is helpful 
in computations to discretise weight into logarithmic bins, such that the weight bracket 
[wi,w i+ i) is the range [wi, (1 + A)ti>j). We adopt a notation: x = log(w/w ), for some 
arbitrary weight wo, and use the function u(x) = Qw<p(w). Thus, integrating u(x) over the 
range [xi,Xi + A), returns the total number of individuals in this size range. 

Some further biological details have to be specified to do the numerical analysis; Table 
Q] summarises the information, and Section 13.11 describes this in more detail. We have 
chosen the parameters not for biological realism but in order to highlight the differences 
between the stochastic jump-growth model and the McKendrick-von Foerster equation. In 
particular we have chosen a smaller predator:prey mass ratio than is realistic. 



metabolic rates being made (see also Law et al. 
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term 


meaning 




value 








Fig [2] 


Fig H 


Fig a 


X 


min wt of phytoplankton 











Xh 

u 


min wt of consumers 


2 


2 


2 


X,i 

CI 


max wt of newborn consumers 


2.1 


2.1 


2.1 




wt at start of senescent death 


5 


7 


8 


X 


max wt of consumers 


7.5 


9 


10 


K 


mass conversion efficiency 


0.2 


0.2 


0.2 


B 


preferred pred:prey mass ratio 


e 1 


e 1 


e 1 


A 


volume searched mass - " 


50 


50 


50 


a 


search volume exponent 


1 


1 


1 


a 


width of feeding kernel 


0.3 


0.35 


variable 


/' 


intrinsic mortality rate 


0.1 


0.1 


0.1 


p 


growth of senescent death 


5 


5 


5 




stochastic realisation 








N 


nnmVipr of nhvtonl^inktnn 


25000 


50000 

{J \J \J \J \J 




No 


initial number of consumers 


2000 


4000 




x 


initial upper bd of consumers 


4 


7 




T *-l 


exponent for fixed spectra 


1.3 


1.3 




A' 


weight bracket for stochastic bins 


0.1 


0.1 






numerical integration 








A 


wt bracket for integration 


0.01 


0.01 


0.01 


St 


time increment for integration 


0.0001 


0.0001 


0.0001 



Table 1: Parameter meanings and values used in computations for figures 



3.1 Model specification for numerics 

The numerical results describe an ecosystem with two types of organism: phytoplankton 
which do not feed on other organisms, and consumers which feed on each other and on 
phytoplankton. In more detail, the full range of body weights [x, x) is subdivided into the 
following regions with different ecological properties. 

• [x, Xf,) is reserved for phytoplankton. These organisms are self-supporting; they do 
not change in mass, and do not form part of the dynamics. Their densities are held 
constant, which is equivalent to assuming that, as soon as they are eaten, they are 
replaced. Such organisms have to be present to provide a supply of food for small 
consumers. 

• [xb, Xd) is a range reserved for renewal of consumers, i.e. a range over which con- 
sumers are born. Renewal is essential: without this, consumers would gradually die 



13 



out. Biological realism requires this range to be distinguished from [x,Xb), because 
newborn consumers may grow in size. When consumers leave this range (by growth 
or by death), they are immediately replaced, which amounts to an assumption of 
perfect density-dependent compensation in the nursery. 

• [xd,Xg) is the range in which consumers experience the standard predation, growth 
and death processes described in Section [2j We include in this range intrinsic mor- 
tality at a per-capita rate fi, which takes into account the fact that organisms can 
die for reasons other than being eaten. 

• [x s ,x) is a range in which the per-capita mortality rate of consumers increases ac- 
cording to the function 



where p scales how fast mortality increases beyond size x s . The purpose of this is to 
ensure that consumers cannot continue to grow indefinitely, in keeping with biological 
constraints on body size. The upper bound x is set such that the density of organisms 
at this size is very close to zero. 

For numerical studies, the predation-rate function k(x, x') needs to be made explicit. In 
keeping with (|T3l) . this function is taken to consist of a volume searched per unit time by 
predators, together with a feeding preference function, which is assumed to have a Gaussian 
shape. In logarithmic variables, the function is: 



where parameters A, a, B remain as defined in Section 12.6} and a measures the range of 
prey sizes likely to be eaten relative to the size of the predator. We have introduced the 
assumption here that predators must be larger than their prey. 

In stochastic realisations, the fixed phytoplankton population was initialised with N p 
individuals taken from an exponential distribution with parameter 7* — 1 over the range 
[x,Xb). The consumer spectrum was initialised with N individuals taken from an expo- 
nential distribution with parameter 7* — 1 over a range [xb,xo). Nq was chosen to make 
the discontinuity between the two spectra small, the upper weight limit being initially Xq 
in the consumers. After the start, consumers dying or growing out of the renewal range 
were replaced with newborn individuals, using the same exponential distribution so that 
the number of consumers in this range would remain constant. We carried out realisations 




/zexp (p(x — x s )) if x > x s 
u otherwise 



(28) 




(29) 
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of the ind i vidua l-based stochastic process (Subsection 12. ip using the Gillespie algorithm 
(iGillespid . Il976l ). Body sizes were aggregated into bins of width A' only for visualisation 
of the size spectra. 



Numerical integrations of the deterministic models were carried out using the explicit 
Euler method, with a bin width A and a time step 5t; consumer spectra were held at 
their initial values in the renewal range. Integrations were initialised with assumptions 
equivalent to those of the corresponding stochastic realisations. For graphical comparison 
with stochastic results, u(x) was scaled such that J u(x, 0)dx was N p and Nq for the 
phytoplankton and consumers respectively, and displayed as the number n(x) = u(x)A ; 
over size intervals A'. 

3.2 Travelling waves 

Figure [2] compares time series from the deterministic jump-growth equation ()12p and from 
the McKendrick-von Foerster equation (|14|) against a realisation of the stochastic process. 
Parameter values are the same for all three time series, and were chosen to contrast the 
two deterministic models, by making the difference between predator and prey body sizes 
relatively small, and by ensuring the steady state would not be an attractor. Initial condi- 
tions were chosen well away from the steady state, to induce large oscillations in the size 
spectra from the start. 

Large sustained waves in density develop over time in all three models. These waves 
move along the size spectra from small to large body size as organisms grow. Peaks of the 
waves are associated with slow growth (prey relatively rare) and low mortality (predators 
relatively rare). As expected, the deterministic jump-growth time series gives a better 
match to the stochastic series than the McKendrick-von Foerster one, in terms of the 
period and shape of the waves (although they are not identical). 

3.3 Variable growth 

The jump-growth model and the McKendrick-von Foerster equation differ in that the 
former describes a process in which organisms, starting at the same weight, develop different 
weights over the course of time. In so doing, the jump-growth model captures an important 
feature of growth: when two organisms of the same weight eat prey items of different 
weights, the two organisms must subsequently have different weights. 

Figure [3] illustrates this feature of the models, using parameter values that highlight the 
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Figure 2: Size spectra expressed as logarithm of numbers logn(x) with logarithm of weights 
x over time t, constructed from (a) the stochastic jump-growth process, (b) the deterministic 
jump-growth equation, (c) the McKendrick-von Foerster equation. Parameter values given in 
TableU 



differences between them. The results show the fate of a set of organisms that all start 
with very similar weights in the range [x^, Xd + A'); the set can be thought of as a cohort 
which grows older without renewal. In the stochastic jump-growth model, organisms were 
tagged individually, and the size distribution of the cohort over time was monitored. In 
the deterministic jump-growth model we assumed a tagged cohort u*(x) at a density low 
enough relative to u(x) for changes in u*(x) to come just from feeding on and being fed 
upon by u(x), without any reciprocal effect of u*(x) on u(x). In the McKendrick-von 
Foerster simulation, differential equations for survival and growth in weight in the cohort 
were solved using the growth and death rates (IT5|) and (TTrjj) respectively, as described in 



Law et al. (120091) . 



The stochastic realisation (Figure [3^l) shows the number of tagged individuals declining 
as time goes on (they are being eaten by larger organisms); it also shows the distribution 
of body weights spreading out. The behaviour of the deterministic jump-growth equation 
matches the stochastic cohort closely (Figure [3h>). However, the McKendrick-von Foerster 
equation (Figure Eb) retains its initial spike-like distribution, because the growth trajectory 
from any size is fixed. 

The average growth trajectories of all three models are close together (Figure [3ji). As 
time goes on and the number of individuals in the stochastic cohort becomes small, fluc- 
tuations in the stochastic growth trajectory can be seen. Also, growth according to the 
McKendrick-von Foerster equation is slightly slower than in the deterministic jump-growth 
equation. However, these differences are small, and it is only when the second moments of 
growth are considered that the spreading in body sizes, missing from the McKendrick-von 
Foerster equation, becomes evident. 

Adding the second-order diffusion term of f|T8|) to the McKendrick-von Foerster equation 
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(a) (b) 




Figure 3: Number n{x) of organisms with log weight x over time t in tagged cohorts embedded 
in size spectra. Cohorts start in a weight range 2.1 < x < 2.2 at t = 0. (a) Stochastic jump- 
growth process; (b) deterministic jump-growth equation; (c) McKendrick-von Foerster equation; 
(d) mean weights over time computed for the cohorts shown in (a), (b), (c), and labelled corre- 
spondingly. Parameter values given in Table [TJ 

f TT^|) would recover the tendency for body size to spread. However, this still leaves out higher 
order terms of the Taylor expansion ( IT8|) which do not necessarily become small unless the 
steady state is an attractor. 



3.4 Dynamical stability 

Figure @] gives examples of the steady states and stability properties of the jump-growth 
and McKendrick-von Foerster models. The breadth of diet a decreases from top to bottom 
in the figure. 

At steady-state, the size spectra have similar shapes in the two models, and diet breadth 
has little effect on them. For the most part the steady states are close to linear under the 
log transformation of both axes. This linearity applies until near x — 8, where the extra 
size-dependent mortality starts to take effect. In the region 2.1 < x < 7 which is close to 
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jump-growth 



McKendrick - von Foerster 




Figure 4: Steady-state size spectra (dashed lines), and transient size spectra (continuous lines) 
after a period of 5 time units has elapsed starting from the same initial function. Column 1 (a, 
c, e) obtained from the deterministic jump-growth equation; column 2 (b, d, f) obtained from 
the McKendrick-von Foerster equation. Diet breadths a: 0.5 (a, b), 0.4 (c, d), 0.3 (e, f); other 
parameters given in Table [TJ Steady state s obtained by New ton-Raphson iteration, which also 



gives the Jacobian matrix at steady state Kess et al I M : numbers given for each graph are 



max(i?e(A)) where A is an eigenvalue of the Jacobian. 
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linear, the slopes are approximately —1.42 in the deterministic jump-growth equation and 
—1.47 in the McKendrick-von Foerster equation, equivalent to exponents 7 = 2.42 and 
7 = 2.47 respectively. These values are close to the value 2.47 predicted from analysis of 
the delta-function version of the feeding preference equation 



Figure H] shows the existence of a bifurcation point at which the system flips from one 
dynamical regime to another as a changes. For large enough a the steady state is an 
attractor, i.e. the Jacobian matrix evaluated at the steady state has max(i?e(A)) < 0: size 
spectra initialised away from the steady state move towards it. For small enough a this 
ceases to be the case, i.e. max(i?e(A)) > 0; instead, the size spectra develop travelling 
waves like those seen in Figure [2j Importantly, the bifurcation point occurs at a smaller 
value of a in the jump-growth equation. This may be because of the lack of spreading in 
body size in the McKendrick-von Foerster equation: such spreading would tend to dampen 
oscillations. A consequence of the difference is that a stability analysis of the M cKendrick- 
von Foerster equation could be misleading; see for instance Law et al. (120091 ). Although 
not shown here, the bifurcation to a travelling wav e can also be ind uced by increasing the 
preferred ratio B of the predator:prey body mass (lLaw et al. i bood i 



4 Discussion 



The starting point for our analysis was a simple, mechanistic, stochastic process in which a 
larger organism feeds on a smaller one, thereby causing the death of the prey and increment 
in its own weight. From the master equation of the process, a macroscopic model for 
the dynamics of size spectra was derived, which we call the deterministi c jump- growth 
mode l. The equation is related to the Smoluchowski coagulation equation f Smoluchowskil . 



19161 ). which describes how the size-distribution of inanimate coagulating particles changes 
over time. However, the jump-growth equation has to deal with special features of living 
organisms, such as their ability to choose the size of their prey, and their inefficiency in 
turning these prey into their own body mass. To cope with the vagaries of the animate 
world, the deterministic jump-growth equation is necessarily more general. 

The expression for the steady-state derived from the deterministic jump-growth equation 
is consistent with the approximate constancy of biomass in logarithmic intervals of body 
mass often observed in marine ecosystems. So the basic empirical regularity evidently 
follows from the bookkeeping of biomass, as it passes through the ecosystem. However, the 
steady state may or may not be an attractor. As one might anticipate from the general 
oscillatory nature of predator-prey systems, another non-equilibrium attractor exists, here 
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comprising waves of abundance that travel from small to large body size. These waves have 
practical as well as theoretical i nterest in view of the large, often unexplained, fluctuations 



in ex ploited marine fish stocks f lHsieh et al 



2006 



Anderson et al. 



2008 



Blanchard et al 



2009 



personal communication). 



The jump-growth model is not the same as the McKendrick-von Foerster equation widely 
used in the study of dynamic size spectra. This is because it allows organisms, starting at 
the same size, to become different through eating prey of different sizes. The McKendrick 



von 



berster equation, with its roots in age distributions ( IMcKendrick 



1926 



von Foerster 



19591 ) does not allow this: organisms which start at the same age must always remain the 

berst er equation has been extended to 



20071 ) , but this was by adding variability 



same age. An age- depende nt McKendrick-von '. 
allow for variable size at age fjGurnev and Veitch , 
to a specified model of growth, the von Bertalanffy equation. Growth of organisms in 
dynamic size spectra comes about in a quite different way, because it emerges directly 
from the action of predators feeding on prey. This is not to suggest that variation in prey 
size is the only cause of variation in predator size; in reality, a variety of extrinsic and 
intrinsic factors are most likely involved. 

Although the deterministic jump-growth model is different from the McKendrick-von 
Foerster equation, the latter can be derived from it using the lowest-order terms in a Taylor 
approximation. The approximation requires that prey size is small relative to that of the 
predator, which will often apply in practice. Thus for many purposes the McKendrick- 
von Foerster equation should work well, notwithstanding the numerical examples used in 
Section [3] (deliberately chosen to contrast the two models). This is with the caveat that 
higher-order terms in the Taylor expansion are not necessarily small when the attractor is 
a travelling-wave rather than a steady state, or when looking at spiky perturbations away 
from the steady state, even if prey are much smaller than their predators. To describe such 
non-equilibrium solutions accurately, the jump-growth model is needed. 

When solving the jump-growth equation numerically, some care is needed in the discreti- 
sation of \ogw. Unlike the McKendrick-von Foerster equation, there is no guarantee that 
feeding will generate non-zero rate terms for growth. If the multiplicative weight brackets 
A are too large relative to prey size, weight increments from feeding do not register, and 
an erroneous solution is obtained. For a Gaussian feeding preference function (12 9|). a rule 
of thumb is that A needs to be of an order Kj (Be 2cr ) to capture properly the rate term 

the order B = 10 2 , a = 0.5 log B and K = 0.1 ar e 



ues o 



1966 



Cohen et al. 



1993 



Jennings and Mackinsonl . 120031 ) 



due to g rowth of organisms. Va 
realistic ( jPaloheimo and Dickid . 
requiring A to be of an order 10 -5 . With marine size spectra encompassing ten orders of 
magnitude, numerical analyses clearly become demanding. A small value of B was used 
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for the illustrations in Section El but it would be much harder to do the computations in 
a more realistic setting. 



Faced with this difficulty, a halfway house would be to use the McKendrick-von Foerster 
equation with the diffusion term from the expansion in f|T8l) . We are not aware of a 
previous derivation of the diffusion term for growth in body size, although diffusion in 
physical space has been c onsidered in the context of the McKendrick-von Foerster equation 
(jOkubo and Levinl . 120011 ) . Nor have we seen the use of a diffusion term in the McKendrick- 
von Foerster equation applied to size spectra, although the effects of introdu cing variability 
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into Gompertz and von Bertalanffy growth models have been described (IBardos 
Gurney and Veitchl . 120071 ). It would be instructive to know how much the McKendrick- 
von Foerster approximation could be improved by introducing this extra term. 

Several further features of real-world ecosystems, not dealt with in this paper, will modify 
our results. First, some feedback between the abundance of phytoplankton and consumers 
is to be expected. Second, perfect compensation in renewal of consumers is unlikely, 
especially when travelling waves affect the abundance of reproducing individuals. Such 
processes generate long, potentially destabilizing, feedback loops. Third, consumers do not 



all start life with the same potential for growt h and reproduction 



of dif ferent species with different life histories (lAndersen and Bever 



hey comprise a number 



2006 



Blanchard et al. 



20091 ). They are born at different sizes, they grow to different sizes, and they allocate 
different propo rtions of their limi ted resources to growth, maintenance and reproduction 
along the way (IMaurv et all 120071 ). Such processes loosen the dynamical coupling between 
a feeding organism and its prey. 

There is much to learn about the intricacies of biology that can stabilize and destabilize 
marine ecosystems. It is important to obtain this knowledge because the biomass in such 
ecosystems is typically of major economic importance, heavily exploited, and with dynam- 
ics that are not well understood. The deterministic jump-growth equation derived here 
should place this programme of research on a more rigorous footing. 
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Appendix A: Derivation of Langevin equation 



Our treatment of the jump-growth model using the van Kampen expansion in Section 12.31 
did not provide a justification for assuming that the fluctuations £ around the solution 
4> of the deterministic equation fl6]) are damped by a factor of Vt 1 ^ 2 . In this appendix we 
derive an approximate stochastic differential equat ion fo r the jump-growth model, adapt- 



ing an approximation procedure used by Gillespie (120001 ) for stochastic models of chemical 
reactions. We will find that the deterministic part of the equation coincides with our deter- 
ministic jump-growth equation §6§ and that the stochastic noise term is indeed suppressed 
by a factor of Vt 1 ^ 2 . 

Because of the stochastic nature of the jump-growth model, the vector of numbers 
"[..., n_i, no, rix, ... ] in each weight bracket introduced in subsection (12. 2 p is described 
by a stochastic process n(t). In a time interval [t, t + r] a number of predation events will 
take place, each of which changes the numbers. This is expressed by the equation 

n(t + t) = n(t) +^ J R ii (n(t),r)i^, (30) 

where the R^ (n, r) are random variables giving the number of predation events taking 
place in the time interval [t, t + r] that involve a predator from weight bracket i and a prey 
from weight bracket j. The are the vectors that give the change in numbers caused 
by such a predation process, as described in Subsection ( 12 .2p . We now will argue that the 
random variables Rij(n(t),r) can be approximated by normally distributed variables. 

The rate of each individual predation event depends on the numbers of individuals 

aij(n) = Vr l kijniTij. (31) 

As the numbers change after each event, the events are unfortunately not independent. 
However, because the numbers change only by ±1 in each event, the change to the rates is 
very small if the numbers are large. So, if we choose the time span r small enough so that 
not too many predation events take place, the rates can be approximated as remaining 
constant throughout the time interval, 

aij{n{t')) « %-(n(t)) W e[t,t + r]. (32) 

In that case the predation events can be treated as independent and therefore the number 
Rij(n(t), t) of event taking place in the time interval follows the Poisson distribution with 
parameter raij{n{t)). 
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Next we assume that the parameter TCLij(n(t)) is either zero or large enough so that 
the Poisson distribution is well approximated by the normal distribution with mean and 
variance both equal to Tdij(n(t)). Again this is easy to justify when the numbers are large 
and provided the feeding kernel fcjj is bounded away from zero. In our case, where the 
feeding kernel contains a Gaussian, we need to neglect the rare events in the tails of the 
Gaussian. 

Note that we are placing two opposing conditions on the size of the time interval r: it 
needs to be both small enough so that the rates don't change much but also large enough so 
that the number of events can be taken to be normally distributed. Such an intermediate 
range for r will exist, provided the numbers of individuals in the weight brackets are large 
enough. In our application, where the overall number of individuals involved is truly huge, 
our approximations will be very good except for very large weights where the density is 
very small and where the approximations will break down. 

Now that we have argued that the Rij are well approximated by normal random variables 
with mean and variance both equal to ra i:) (n(t)), we express them as 



Rij(n(t), t) = a ij (n(t))T + yo^(n(i))r ry (33) 

where the are normal random variables with mean and variance 1. Substituting this 
into (130|) . rearranging terms and dividing by r gives 

n(t + r)-n(t) = ^ a .^ n ^ u .. + ^ ^ a .. {n^u^r' 1 ' 2 ^. (34) 



T 

ij ij 



We now approximate this equation, which is valid for small but finite r, by the stochastic 
differential equation obtained by taking the limit r — >■ 0, 

~~rr^ = "■,{'»■{' + \Z a ij{ n (t))»ijVij{t), (35) 



where are independent white noise process e s. T 



Langevin equation, see for example Ivan Kampenl ( 119921 ). 



lis type of equation is known clS cL 



Substituting the explicit expressions ( 13Tj) for the rates into the Langevin equation ( 13"5]) 
gives 

— — =f2 ^ ^ ( kijUjUj kjiUjTii + k m jn m Tij^ (36) 
j 

+ / ^ ^ ^ — "\f kij Tii TLj Tjij — ■sj kjiTijTliTjji + k m j Tl m 71 j 7] m j ^ . 
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When we write the equation in terms of the population densities $j = rii we see that 
the fluctuation terms are suppressed by a factor of Vr 1 ^ 2 . 



j 

For large system size Q the fluctuation terms can be neglected and we end up with our 
equation (jBJ). 
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